Combined computational modeling and experimental study of the biomechanical mechanisms of platelet-driven contraction of fibrin clots

While blood clot formation has been relatively well studied, little is known about the mechanisms underlying the subsequent structural and mechanical clot remodeling called contraction or retraction. Impairment of the clot contraction process is associated with both life-threatening bleeding and thrombotic conditions, such as ischemic stroke, venous thromboembolism, and others. Recently, blood clot contraction was observed to be hindered in patients with COVID-19. A three-dimensional multiscale computational model is developed and used to quantify biomechanical mechanisms of the kinetics of clot contraction driven by platelet-fibrin pulling interactions. These results provide important biological insights into contraction of platelet filopodia, the mechanically active thin protrusions of the plasma membrane, described previously as performing mostly a sensory function. The biomechanical mechanisms and modeling approach described can potentially apply to studying other systems in which cells are embedded in a filamentous network and exert forces on the extracellular matrix modulated by the substrate stiffness.


Supplementary Note 1 -Model description
In this section we provide some additional details about several aspects of our multiscale computational model. Also, Supplementary Table 1 includes a list of processes and space and time scales represented in the multiscale model.
The model requires the use of several computational techniques to avoid unreasonable runtimes. First, all independent operations are performed in parallel using the Thrust parallel algorithms library to optimize block and thread counts. Force calculations for individual platelet filopodia, fiber nodes, and platelet-fiber interactions are performed on Nvidia GPUs with synchronization before and after solving the respective ordinary differential equations in time. Second, since all filopodia and fibers must interact, two 3-dimensional cell-grids of different resolutions are utilized to build nearest-neighbor maps.

Supplementary Table 1: Clot contraction processes and scales Processes
Spatial and temporal scales Microscale: Fibrin accumulation near platelets Spatial: from a few to tens of µm (depending on platelets cluster size) Temporal: from a few minutes to a few hours (30 minutes in the model simulations) Microscale: Fibrin-platelet interactions Spatial: from less than 1 µm to a few µm Timescale: from a few seconds to a few minutes Mesoscale: Platelet clustering Spatial: from few to tens of µm (depending on platelet cluster size) Temporal: from a few minutes to a few hours (30 minutes in the model simulations) Macroscale: Fibrin network contraction Spatial: from tens of microns to mm (47.6 µm maximum clot diameter in our simulations) Temporal: from a few minutes to an hour (30 minutes in the model simulations) Macroscale: Contraction of a whole clot Spatial: from tens of microns to mm (47.6 µm maximum clot diameter in the model simulations) Temporal: from a few minutes to an hour (30 minutes in the model simulations) The grid for the network has a lattice size given by 10 ⋅ 2 !"# = 1 while the grid for the platelets has a lattice size of 1.5 ⋅ !"$ = 3.3 . This grid method reduces the cost of element interaction from ( % ) to ( ), resulting in scalability and nearly linear runtime increase with respect to element count. In our computational model, both main components (fibrin network mechanics and platelet mechanics) act on mechanical properties and generate mechanical forces.
At each time step the neighbor list for each fibrin (node or sub-node) and for platelets is built using the grids described above. Only the nodes in the same neighborhood are considered when calculating interactions. This is done to reduce the computational cost while maintaining high computational accuracy. The forces acting on each node (fibrin node or sub-node or platelet) are then calculated. In particular, at each time step the strains of all fiber sub-segments are calculated and used to determine the elastic forces, between fibrin nodes and sub-nodes, and to determine the pulling force of each currently pulling filopodia as their force depends on the strain of the sub-segments connected to the node (or sub-node) that is currently being pulled. Then a forward Euler method is used to solve the system of differential equations for fibrin fibers and platelets nodes. Finally, the process is repeated for the next time step until the prescribed computational time is exhausted.
Supplementary Note 1.1 -Individual platelet and fibrin fiber sub-models. Platelets and fibrin dynamics are described by Langevin equations as follows. Motion of the location of each platelet center of mass is represented by a solution of the stochastic ordinary differential equation where &," is the deterministic force acting on the -th platelet, &" is the viscous damping force, and " ( is the Brownian force due to thermal fluctuations 1 . We consider the system to be in an overdamped regime as in Britton et al. 2 , due to the small Reynolds number. The inertial term is therefore neglected, and the system (1) is discretized in time using a Forward Euler scheme: where the superscripts n and n+1 refer to the vector quantitates at time steps n and n+1 for the -th node and = 0.001 is the time step. & represents the drag coefficient for a platelet with fixed volume and it is calculated using the Stokes estimation, & = 6 & 3 . The blood viscosity is = 0.004 . and the platelet radius is & = 1. 13 , which results in the value & = 85.2 . . /+ . Similarly, dynamics of fiber nodes is determined by solution of a Langevin equation where !,0 is the deterministic force acting on the -th fiber node, and ! is the drag coefficient for a fiber node with mass !,0 . The system is discretized and solved in time using the same scheme: Assuming fibers have a diameter of 100 nm, the value of ! was calculated using & = 0.05 to be ! = 3.77 . . /+ .
The Brownian force in three dimensions is calculated using the Einstein relation between the diffusion coefficient of a medium and the temperature 4 . Because the mean squared displacement of a particle following Brownian dynamics is calculated as . over a single time step, the random force acting on a fiber node is " , where " is sampled from a standard normal distribution.
Random forces are calculated similarly to the case of platelets with the drag term, & . The deterministic forces, &," and !,0 , acting on each platelet's center of mass or fiber node, are separated into components as follows. For each platelet, the forces are generated by (1) pulling forces of filopodia, (2) adhesive forces and (3) volume exclusion forces. The total force is calculated using the negative gradient of the energy, " = − " 453 , where " 453 is the energy associated with the -th platelet separated into three components: " 453 = " 675 + " 89: + " ;<= (5) with " 675 , " 89: , and " ;<= representing the filopodia, adhesion, and volume exclusion energies determined as follows: where = " − 0 is the distance between the -th platelet and the -th platelet/fiber, = 2 A & /2 +/1 , and = @ % ⋅ 1 / = 16 ⋅ . Note that in the main text that we refer to | − 675 | = &$>-?$?-, − 89: = @ , and − ;<= = 5C . Since we did not have specific data to calibrate the Lennard-Jones forces, we have chosen parameter values based on computational stability. The minimum of the Lennard-Jones potential (see formula 8) is achieved when = . Our choice of parameter values places the minima at = = 2 A & /2 1 < 2 A & assuring that the minimum of Lennard-Jones energy is achieved slightly inside the incompressible bodies of two neighboring platelets (since platelet bodies can partially deform when they aggregate). The choice of = @ % was determined by the assumption of Lennard-Jones forces being stronger than the adhesive force @ and computations showed stable results using = @ % . Similarly, the force acting on each fiber is calculated using the same formula (5) but with the force term separated into (i) forces applied externally to the network by platelets which generate (ii) individual fiber forces and (iii) bending forces. Each individual fiber is represented by a nonlinear Worm-Like-Chain (WLC) spring and by bending springs introduced at each node triplet as in Britton et al 2 . Thus, the total deterministic energy of the -th fiber is as follows: where 0 ;<3 , 0 D5= , and 0 (;E9 represent the fiber, WLC and bending energies, respectively: where # is the Boltzmann constant, is the absolute temperature, and are the scaled persistence and contour lengths. The normalized strain, ̄ is calculated between the -th and -th nodes 5 . Each spring is represented using for the number of fibrin monomers in parallel, and the normalized contour length C representing the number of fibrin monomers in series 5 . The energy of WLC springs acting on node is calculated by summing up impacts from neighboring nodes. The energy on the node is calculated by summing over all triplets, ( ), surrounding it. Bending energy is equal to the bending constant, K , multiplied by the change between the current angle $02 and the preferred angle @ between nodes , , and (see also Figure 3 in the main text).
Formation of cohesion bonds leading to fiber crisscrossing in our model, instead of the volume exclusion principle, prevents fibrin fibers from passing through each other. Specifically, at each time step, distances between initially disconnected fiber nodes are compared with a threshold. If the distance between nodes of different fibers is less than the fiber diameter of 0.1 μm, a flexible link is established between the nodes to represent a fiber-fiber cohesive bond (See also Britton et al. 2 .) Moreover, as the platelet's filopodia pull fibers towards the platelet body, fibrin accumulates on the surface of the platelet. Volume exclusion between platelet and fibrin fibers is implicitly represented in the model as follows. Once the surface of the platelet gets saturated by fibrin as a result of filopodia pulling on fibrin fibers, any nodes representing additional fibrin fibers near platelet surface will be repelled by the WLC forces. This repulsion eventually balances with the filopods' pulling forces, leaving the clot at rest. This accumulation of fibrin on the platelets surface can be seen in Figure 7.

Supplementary Note 1.2 -Platelet-platelet interactions. Two platelets can interact if the distance between their centers of mass
the respective regions 1 of the two platelets overlap but do not overlap with region 2 or 3 (see Figure 4 in the main text). If the filopod attaches to another platelet, when 2 & < && < 2( & + !"$ ), the force & ( ) = @ /2, the minimum force exerted by a filopod (see formula (5) in the main text). This model assumption was chosen so that if two platelets are pulling on each other the total force exerted on each platelet will be @ (since @ /2 will come from the reaction force of the first platelet pulling on the second and another @ /2 would be exerted by the other platelet). If the distance between two platelets is such that A ⋅ 2 & < && < 2 & , the two platelets will be in each other's adhesive zones. Note that in our simulations 45 platelets are initially uniformly distributed in space at an average distance of 25 microns from each other.

Supplementary Note 1.3 -Variations in the number of filopodia per platelet.
The number of filopodia for each platelet in specific clot contraction simulations is extracted from a log-normal distribution based on the experimentally observed numbers of filopodia per platelet ( = 1.85, = 0.27, see Figure 5 and the Results section in the main text) and two more distributions with a lower number of filopodia (lognormal with = 0.9 and = 0.19) and a higher number of filopodia (log-normal distribution with = 3.7 and = 0.38. The distribution of platelet filopodia count with larger corresponds to a higher degree of platelet activation, while the distribution with smaller values of corresponds to a lower degree of platelet activation. When initializing the number of filopodia for a specific platelet in our simulation a random number is generated using one of the distributions mentioned above (e.g. normal, low, or high number of filopodia) and if the number is between 2 and 13 (minimum and maximum number of filopodia per platelet) the integer part of the number is used to determine the number of filopodia for that platelet in our simulation.

Supplementary Note 1.4 -Parallelization and model computational cost.
The GPU computational cost remains constant with the increase of the number of fibrin fiber nodes. This is because each node is assigned with a GPU thread to update its position. Since the GPU computing uses the shared memory model as opposed to the CPU parallelization in which the distributed memory model is used, the parallel communication time of the GPU computing in theory is zero. This is a significant advantage over a CPU implementation in which MPI is typically used for message passing (data transfer). The chief computational bottleneck is the GPU memory size. We note that the limitation of the current GPU implementation is the size of the physical memory of the GPU card. If the fiber network data must be split to store on multiple GPU cards, then additional parallel communication time overhead for communication between different GPU cards is introduced. However, the number of nodes in this model simulation study is on the order of thousands or tens of thousands of nodes and can be stored on a single GPU card.
Moreover, the limited number of CPUs endemic in CPU parallelization would mean that, rather than assigning a single node to each process thread for the relatively simple calculations of forward Euler forces (as is the case in GPU parallelization), optimization of CPU-parallelized code would involve assigning nodes to groups and processing each group of nodes as a CPU thread. The assignment of these groups would not necessarily be straightforward, as it would be platform-specific in its implementation based on the number of available cores and decrease code readability without a clear computational advantage since, in effect, a more granular (node-level) parallelization is natural to this model.

Supplementary Note 2 -Initiation of structure of a fibrin network with platelets
The Supplementary Note 2.1 -Platelet and fibrin nodes placement. First, a spherical domain of radius, , is set. Next, a fiber density, ! , and platelet density, & , are chosen. These parameters are used to calculate the total number of fibers and platelets. The platelet effective radius, 4 , is chosen and a fixed fiber density parameter is chosen which represents the percentage of the total number of fibers from the distribution 0 < < 1 that will be created. In the algorithm, we set = 0.975, resulting in less than 2.5% error in the node degree and fiber density. Platelets are uniformly distributed inside the domain using Latin Hypercube Sampling (LHS) 6,7 . If a platelet, 0 , is placed outside the spherical volume in the LHS sampling procedure or if the platelet is far from all fiber nodes, new uniformly distributed coordinates are chosen. The new coordinates are regenerated until " ( ( 0 , " )) < 4 , and ","L0 ( ( 0 , " )) > 4 ,where the minimum is taken over all fiber nodes, " , and platelets, " , respectively. These conditions ensure that all platelets are near some fiber point and far from all other platelet points. Once platelets are placed in the domain, the 4 fiber points are placed so that each point is randomly placed near a platelet within half the distance of the platelet radius, 4 , and the 6 fiber points are randomly placed throughout the spherical domain.
Supplementary Note 2.2 -Degree distribution error reduction. Next, fiber nodes are given a preferred degree, " & , by sampling the distribution, . When fibers are generated between nodes, a node is considered to be in need of a neighbor if its current degree, " M , is such that " & > " M . Fibers are created " M according to the following algorithm. A node, " , with " & > " M is chosen at random to start a fiber. An end node must then be chosen to complete the fiber. All points, 0 , with " = ( " , 0 ) < are weighted with weight ( " ). Next, a single sample is drawn from . The node 0 satisfying | ( " ) − | is chosen to end the fiber if the fiber between " and 0 has not yet been established. If the fiber has already been established, we set ( " ) = 0 and choose a new sample * from , resulting in a new end point for the fiber beginning at " . If no node can be found, a new random starting node " is chosen again from the nodes such that " & > " M . This process of creating fibers continues until the desired density of fibers is reached. Choosing fibers in this manner serves to approximate the target distributions for experimental fiber length and node degree.

Supplementary Note 2.3 -Length distribution error reduction.
To increase the accuracy between the current fiber lengths and{ K }, nodes are perturbed randomly within the domain. A random move serves to slightly change the associated fiber lengths. If a point from 6 is chosen, any move is allowed within the fixed domain. If a point, " , from 4 is chosen, the random perturbation must satisfy ( " , ) < 4 , where is the platelet associated to " . To determine fiber length error, a sample is drawn representing each fiber from the experimental distribution. Samples are binned with width (=0.1) and errors are determined by comparing the frequency in each bin. Using this method, a node perturbation is accepted if the new fiber lengths lower the fiber length error.

Supplementary Note 3.1 -Calibration of an individual fiber sub-model. Fibrin fiber elastic properties
were calibrated using the general WLC modeling approach. The WLC model has been previously shown to provide analytical and theoretical explanations for single fibrin fiber mechanical behavior 8 . Here, we use the experimental stress-strain curves for a single fibrin fiber obtained using atomic force microscopy (AFM) 5 to calibrate each simulated fiber response to stretching. Two parameters, persistence length and contour length multiplier, in the WLC model were determined by using a recursive least-squares algorithm 9 . The experimental data and resulting fitting curve for a single fiber are shown in Supplementary Figure 1.

Supplementary Figure 1.
Single fiber force-strain response calibration. Black line shows the model fit to experimental data for single fiber (blue dots) forcestrain response using an atomic force microscope AFM 5 . Table 2 lists parameter values that are calculated based on the recursive least-squares optimization under the assumption that the temperature is fixed at 300 Kelvin and 200-1100 monomers would fit in the circular cross-sectional area of a fiber 5 . To incorporate fiber bending, bending springs are added between sub-nodes for all fibers (see Figure  3B in the main text). A bending spring is placed at every node triplet (see Figure 3B in the main text). Deviation from an equilibrium angle results in generation of a restoring force that is calculated with the bending spring potential (Eqn. 11). The bending spring constant, K , is calculated using the relation K = , where is the Young's modulus based on the linear response of individual fibers, and is the moment of inertia for a circular shape. Assuming all fibers maintain the same diameter of 100 nm, we can calculate Based on data from Lam et al. 11 , the parameters in the equation (12) are obtained through least square fitting as follows. First, we used the data in Lam's paper Figure 2C and a minimum force of 4 at @ (no strain) to fit , , , @ using least square fitting with &$>-?$?-( , ) = @ + >2 #2*M that is obtained from formula (12) for platelet force at =∞, when the cantilever reaches equilibrium (see also 11 ). This resulted in the values = 61, = 3.6, = 130, and @ = 2.2 (see also Supplementary Figure 2 and Supplementary  Table 3). Second, we use the data in supplementary material of Lam's paper figure 1A with = 18 / (chosen because the typical range of stiffness in our simulations is [2.71, 20.81] / ) and the values of , , , @ previously determined to find and obtained using least square fitting for formula (12). This resulted in the values = 1.4 and = 640 (see also Supplementary Figure 2). However, equation (12) cannot be used as is since Lam's data refers to a platelet attached on two sides to a cantilever and a stiff surface. To adapt it for use with filopodia we proposed the following modifications with & ( ) representing the number of filopodia pulling at time . Equations (13) and (14) both assume that the measured data of platelet force from Lam et al. 11 comprise forces on two sides of the platelet, and that only a single filopod will act on each side (hence the division by 2). However, equation (13) also assumes instantaneous force loading time for the filopodia. Similarly, equation (15) also assumes instantaneous loading time but does not divide the force by 2. Finally, equation (16) assumes that the total force of the platelet of equation (12) is instantaneously loaded and equally divided between active platelets at time .  Figure 6). The transition from the full-platelet force expression (12) to the single-filopod expression (13) is also mechanistically reasonable. In the experimental settings of 11 , the force on a platelet has been measured between the two broad surfaces of a cantilever and the substrate, while the filopodia can have only one point of adhesion (i.e., one filopod to one fiber). Therefore, the force measured in these experiments may be the result of the platelet exerting force with two filopodia -one attached to the surface and another to the cantilever. This suggests that the force for one filopod should be expressed as half of what is given in formula (12), if indeed (12) was the magnitude of force measured from two filopodia.

Supplementary
Furthermore, during active contraction, platelet filopodia contract much faster than the platelets in the atomic force microscopy experiments 11 . Therefore, the model assumes a fast loading time by setting the time t equal to infinity in equation (12) and dividing the force by two because the filopodia has only one focal point for the exerted force instead of two like in 11 .  Figure 7C is fitted with a polynomial function and then the most relevant peaks (non-negligible maxima) of the fitted polynomials are selected as the beginning and ending time of the kinetic clot contraction phases. The exact numbers are reported in Table 1 of the main text together with the average contraction rates in the phases calculated as the mean change in fiber density during each phase.

Supplementary Note 4.2 -Impact of the number of platelet filopodia on the clot contraction kinetics.
The change in fiber density over time from Figure 8C is fitted with a polynomial function and then the most relevant peaks (non-negligible maxima) of the fitted polynomials are selected as the beginning and ending time of the kinetic clot contraction phases. The exact numbers are reported in Table 1 of the main text together with the average contraction rates in the phases calculated as the mean change in fiber density during each phase.

Supplementary Note 4.3 -Clot contraction kinetics.
Platelets have been observed to sense and dynamically respond to variations in the clot microenvironment, such as alterations in substrate stiffness and locally applied forces [11][12][13] . While individual cellular reactions may vary, the general consensus reveals a correlation between increased substrate stiffness and increased expression of the active integrin 3, platelet spreading, and contractile force 12 . In this section, we present several different simulated responses of platelets due to changes in the stiffness of the extracellular matrix and compare the contraction metrics obtained, with a particular focus on the densification rate (see Supplementary Figure 6C) and the contraction phases observed in experiments. In particular, we tested variations of the filopod force obtained in Supplemental Note 3.2. We compared simulations with filopodia forces given by equation (13) and the other force types obtained from the general equation (12) (see Supplemental Note 3.2). These forces were obtained for different hypothesized force loading rates (i.e., how long it takes a platelet from the beginning of contraction to achieve full pulling force) or concerning the division of the platelet force among active filopodia (i.e., are the filopodia independent or dependent on each other).  (Supplementary Figure 6). Moreover, for these simulations, the rate of fiber densification revealed several phases that qualitatively agree with the densification phases observed in experiments 14 . In particular, in these two simulations, it is possible to observe a growth in the densification rate (phase P1), a decrease (phase P2), and an asymptotic limit (phase P3). These simulations revealed that even similar force profiles (based on equation (12)) can produce significantly different results. Also, we observed that the force described by equation (13)  The speed of platelets during clot contraction was evaluated. These simulation results allow us to determine if the platelet positions are tending toward equilibrium (i.e., speed decreases) as the contraction progresses (i.e., higher simulation time). This relates to the formation of clusters, in that as platelets slow down, it takes longer for them to approach each other and for clusters, and therefore platelets clusters form less at later simulation times. As construction progresses and elastic forces in the fibrin network reach an equilibrium with the platelets pulling forces, the platelets start moving less and the clot approaches equilibrium, as can also be deduced by the slowdown in contraction rates observed in Supplementary Figure 4. Note that as our system is in an overdamped regime a low velocity corresponds to approaching an equilibrium of forces. Supplementary Figure 8 shows the probability density of the speed of platelet relocation during clot contraction at 10, 20, 30, and 40 minutes. As the clot contracts, the fraction of slower moving platelets increases and, as a consequence, platelets form clusters more slowly as clot contraction progresses.  Figure 9 shows the angle distribution in the x-y plane of the platelet clusters in our simulated clot at t=10, 20, 30, and 40 minutes of clot contraction. This angle is calculated with respect to the x-axis in our simulation coordinate axis and its relative angle with respect to the projection on the x-y plane of the line connecting the center of our clot and the center of a cluster of platelets. We can observe that two-platelet clusters are equally distributed with respect to their angle in the spherical clot. Moreover, bigger size clusters, while less in number, do not show any marked asymmetric distribution (see also Supplementary Figure 7). Other 2D projections (x-z and y-z plane) show similar pictures. Supplementary Figure 9 demonstrates that no radial asymmetry forms in the spatial distribution of platelets clusters in our model. Figure 9. Angular distribution of platelet clusters at various times during clot contraction. The angular probability density function of platelet clusters with 2 or more platelets in the x-y plane is shown by a pie chart. Different colors correspond to different time moments of clot contraction in simulations. The black bars represent the error as one standard deviation from the mean indicated with the corresponding-colored histogram bar. Figure 10 shows the evolution of strain in time across multiple simulated clots. Looking at the mean strain per fibrin segment will show only the mean and standard deviation of up to strain 0.1 at approximately 40 minutes. It may be important to note that in Figure 7, we observe fiber strain of up to 1.8 localized into bundles between platelet clusters. Additional simulations were run wherein we reduced the force magnitude that filopodia exert on fibrin to one-tenth of experimentally calibrated values to replicate the effects of a non-muscle myosin IIA inhibitor treatment. Our quantifications of this effect are shown in Supplementary Figure 11.

Supplementary Note 4.5 -Additional analysis of fibrin strain. Supplementary
Supplementary Figure 10. The means and standard deviations of mean-strain-per-node for different types of simulated clots are shown. Each dotted line and ribbon shows the average value and standard deviation of mean-strain-per-filopodia across ten simulated clots. Simulations with experimentally observed intermediate filopodia numbers are shown by the green line, while blue dashed lines and red dotted lines represent average strain with the higher and lower filopodia numbers, respectively.